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A recently developed Quantum Monte Carlo algorithm based on the stochastic evolution of 
Hartree-Fock states has been applied to compute the static correlation functions of a one-dimensional 
model of attractively interacting two component fermions. The numerical results have been exten- 
sively compared to existing approximate approaches. The crossover to a condensate of pairs can be 
identified as the first-order pair coherence extending throughout the whole size of the system. The 
possibility of revealing the onset of the transition with other observables such as the density-density 
correlations or the second-order momentum space correlations is discussed. 
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I. INTRODUCTION 



The recent developments in the cooling and trapping techniques of neutral atoms have opened the way to the realization 
of fermionic atomic samples at temperatures well below the degeneracy temperature pj. This suggests that atomic 
gases are ideal candidates for the study of the physics of degenerate many-fermion systems. With respect to solid 
state ones, atomic systems offer in fact a better isolation from external disturbances such as material defects, a better 
knowledge of the microscopic details of the systems, as well as a wider range of tunability of the parameters, in 
particular the interparticle interactions. By tuning the external magnetic field around a Feshbach resonance, the 
atom-atom scattering length a can be varied from kpa = —00 to +00 (kp being the Fermi momentum) opening the 
way towards a comprehensive study of the pairing transition both in the regime a > in which a Bose-condensate 
(BEC) of tightly-bound molecules is present, and in the regime a < (BCS) in which a condensate of Cooper pairs 
is formed. Diatomic molecules have been created and observed by several experimental groups 0- Bose- Einstein 
condensation of tightly bound diatomic molecules has been recently reported |3| • The crossover region between BEC 
and BCS is currently under experimental investigation 4] and first evidences of pairing in the crossover region have 
been reported in 

From the theoretical point of view, a large effort is currently made to establish the main features of the pairing for 
high values of the scattering length kp\a\ 3> 1, regime in which the atomic gas shows strong correlations [aSU- I n 
particular, the dependance of the transition temperature on the interaction strength in this crossover region is still 
an open problem. 

The present paper reports a numerical study of the condensation of pairs in a regime of relatively strong interactions, 
so to characterize the consequences of the transition on the different observables of the system and identify specific 
features which may represent unambiguous signatures of the onset of condensation of pairs. 

The calculations have been performed by applying the quantum Monte Carlo (QMC) method developed in [|| to 
a one-dimensional lattice model of fermions with attractive on-site interactions. A short description of the model 
under examination is given in sec. [HI while the numerical algorithm used for the calculations is presented in sec. IIIII 
Numerical results are presented in sec . 11 VI and then extensively compared to the predictions of a perturbative expansion 
in the interaction coupling constant (sec. [Vj), and of existing approximat e ap proaches (sec lVI|) . such as the BCS 
theory 111 EH, two versions of the random phase approximation (RPA) [T2, Il3 as we U as the Nozieres Schmitt-Rink 
theory |g. 

Several among the most relevant correlation functions of the Fermi gas have been considered here, in partic- 
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ular the opposite-spin density-density correlation function (pi(x) p-j- (0)) , the first-order pair coherence function 
(^|(a;) ^|(0) ^^(O)) and the second-order momentum space correlation function (nk-\h-ki)- The density- 

density correlation function has been already the object of several papers studying the experimental signatures of the 
BCS transition in atomic Fermi systems, e.g. [l4j|. while the first-order pair coherence function is the counterpart, in 
a non-symmetry-breaking approach, of the order parameter of the phase transition in a Landau-Ginzburg theory |lCj . 



II. THE PHYSICAL SYSTEM 



A one-dimensional low energy two-component Fermi gas can be modeled by the Hamiltonian: 

k,a x 

The spatial coordinate x runs on a discrete lattice of Af points with periodic boundary conditions; L is the total 
length of the quantization box and dx = L/Af is the length of the unit cell of the lattice. The spin index runs over 
the two cr =1,1 spin states. The system is taken as spatially homogeneous, m is the atomic mass, and interactions 
are modeled by a two-body discrete delta potential with a coupling constant g . The field operators ^ a (x) satisfy 
the usual fermionic anticommutation relations {'^> tT (x), (x 1 )} = 5 a . a i 8 x ^ X '/dx and can be expanded on plane waves 
according to ^ a (x) = ^2 k d^e 1 ^ / \f~L with k restricted to the first Brillouin zone of the reciprocal lattice. In order for 
the discrete model to correctly reproduce the underlying continuous field theory, the grid spacing dx must be smaller 
than all the relevant length scales of the system, e.g. the thermal wavelength and the mean interparticle spacing. 
In the present one-dimensional case, the relation between the coupling constant on the lattice and the physical ID 
coupling constant guj is: 

90 = 9iD [I + ^ h2 ) , (2) 

which, in the limit dx <§; tt 2 Y?jmg\u reduces to the expected one go = giD ITEL Hfij. This condition is satisfied 
in the Monte Carlo simulations presented in this paper. We also note that two particles interacting in free space 
with a attractive delta potential in ID have a bound state of energy —mg 2 D /4h 2 . In the numerical examples of this 
paper, the Fermi energy is much larger than this binding energy so that we are not investigating the condensation of 
preformed pairs but rather a BCS regime. 



III. THE QUANTUM MONTE CARLO SCHEME 



We assume the gas to be at thermal equilibrium at a temperature T in the canonical ensemble, so that the unnormalized 
density operator /O cq (/3) = e~^ n with (3 = 1/ksT. From textbook statistical physics, we know that such a density 
operator can be obtained by means of an imaginary-time evolution: 

^^ = -lin Pc(1 (T) + p C(1 (T)H] (3) 

during a "time" interval r = — > (3 starting from the initial state corresponding to the infinite temperature case 
where p C q(r = 0) = ljv, 1 being the identity matrix in the A^-body Hilbert space. 

As it has been recently shown in 0, the exact solution of the imaginary-time evolution © can be written as a 
statistical average of Hartree-Fock dyadics of the form: 

a = |^...^)(^...^|. (4) 

For 01 = 1,2, (t>f ] (j = 1 ... TV) are Hartree-Fock orbitals for the iV fermions, in the sense that: 

|#>...^> = at .at,, | 0)j (5) 
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the creation operator corresponding to the wavefunction <j)(x, a) being defined as: 

4 = £ dxcj>(x,a)¥ a {x). 



(6) 



For the model Hamiltonian (JTJ , the imaginary-time evolution of each of the orbitals <pj can be reformulated in terms 
of Ito stochastic differential equations of the form: 



2m 
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(x, -a)\^ (x, a) - <j>f> (x, ~a) 4>f ] (x, a) ^ (x, -a) 



2 t-J ^ |U( Q )||2 



j x',a' 



^* (x>, a') (x', -a') <pf (x>, -a') ^ (*'. "0+ 



-^ ) *(^ CT 0^*(^,-^)^ a) (^'^')^ a) (^-^)] ^W)} (7) 

where P represents the momentum operator on the grid and the norm ||(/>|| is defined as \\4>\\ 2 = ^2 xrT dx \4>(x, cr) | 2 . The 
deterministic part is simply the mean-field Hartree-Fock equation in imaginary time, while the correlation functions 



of the zero-mean noise dB^ are given by: 

50 



dB < f ) {x,<T)dB { "' ) {x',<j') 



The projector Q "} x CT \ projects orthogonally to the subspace spanned by the wavefunctions (j)^ a> (x, a). A possible 
noise with the required correlation function (JSJ) is: 



fj T n {a) o (q/) 

2dx 



f>i(x, cr) <f)f \x' , cr') 5 a ,a> S a ,-a> S x , x < 



(«), 
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dB^ixA) 



go 
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o e a >*(x) 
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with (x) independent zero-mean Gaussian noises with £( a )(x) ^ a '^(x') = 0, £( a )*(a;) £( Q ')(a;') = 5 x ,x>^a,a>- It can 
be proven that this set of stochastic differential equations reproduces, in the average over the noise, the exact 
evolution of the Hartree-Fock dyadic a during dr: 



da = -y [«cr + <7«]. 



(10) 



The initial state In can be written as a functional integral over all possible sets of orthonormal wavefunctions 
(j = l...N): 



(11) 



This writing of the identity operator can be used as a starting point for an exact simulation of the fermionic many- 
body problem. To this purpose, we have to numerically solve the stochastic differential equations Q for imaginary 
times going from r = to r = (3. This is done by splitting the imaginary-time interval into a large enough number M. 
of time steps; Q a \x) is the noise terms at the time-step j (j = 1 . . . M.) on the site x. The expectation values of any 
observable at temperature T is then obtained as an average over all the possible values of the initial wavefunctions 
(jjf^ and the elementary noise terms £j a \x). 

For example, the partition function Tr[p] is obtained as: 



TrMH^-.^W...^ 



(12) 



or, equivalently, as the determinant Det[M] of the matrix M whose entries are Mij — (<j)p \4>^)- In the following, 
we shall be mainly interested in the one- and two-body correlation functions of the gas. By making use of the Jacobi 
theorem [l7j |. these can be usefully rewritten in the following compact forms: 



*t (*)**<(*')) =Det[M] Y, (M-\<tf*(x,a)$\ X >,i/) 



(13) 
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and 



Det[M] ^2 Dct 

ijkl 



^*(x, a) <f>f*{x>, a') ^ (//) ^ (*"', a'") (14) 



In a practical simulation, the averages are performed by means of Monte Carlo techniques. A description of the details 
of the numerical algorithm used is given in Appendix E] 



IV. MONTE CARLO RESULTS FOR THE CORRELATION FUNCTIONS 



A Monte Carlo code based on the stochastic approach described in the previous section has been used to numerically 
compute the expectation values of some one- and two-body correlation functions for a one-dimensional Fermi gas with 
attractive binary interactions as described by the Hamiltonian with go < 0. The results of analogous calculations 
performed with a very similar Monte Carlo algorithm have been reported recently in 0] . Other Quantum Monte Carlo 
schemes have also been applied to the numerical study of the fermionic Hubbard model with attractive interactions 
at finite temperature. In particular, the determinantal QMC algorithm has been used to study the correlation 
functions in 2D [2(j and the transition temperature to a pair condensate state in 2D |2ll and in 3D [2^ . 

For our simulations, a lattice of Af = 16 points was taken, with a total number of N — 12 atoms. A number M. 
of imaginary-time steps comprised between 400 and 1000 has been used. As already mentioned, the ensemble in 
which observables are calculated is the canonical one; note that the number of particles in each of the spin state can 
fluctuate, only the total number of particles is fixed. As the two spin components are equivalent, the mean densities 
in each of the spin components are equal: 

N 

Pr = Pi= 2L' ^ 

The state of the gas in the absence of interactions and at T = is depicted in Fig. ^ m a given spin component, 
the 5 lowest-lying single particle energy levels are totally filled, whereas the two degenerate states of wavevectors 
k + = kp — 6tt/L and fc_ = — kp are half-filled. More precisely, 10 atoms are frozen in the states of \k\ < kp, and the 
two remaining atoms are distributed among the 4 degenerate states, | f or J., ±kp), which can be done in 6 different 
ways. In presence of attractive interactions, this degeneracy will obviously be lifted and the configurations with one 
atom | and one atom J, with opposite momenta in the degenerate multiplicity are favorable to the formation of a 
Cooper pair. 



A. One body correlation functions 



The simplest observable to compute is the one-body correlation function in a single spin state a (normalized to the 
density p a ): 

gi 1 }(x) = -(^(x)^ a (0))- (16) 

Per 

The Monte Carlo prediction is plotted in figElfor different values of the temperature: as expected, this correlation 
function is short-ranged, coherence extending only on a length of the order of the Fermi length lp = 1/kp for 
T < Tp . This correlation function is indeed the Fourier transform of the momentum distribution of the gas. As the 
interactions affect the momentum distribution only in a thin region around the Fermi surface (the Fermi points in our 
one-dimensional geometry), they do not significantly modify its shape as compared to the ideal Fermi distribution. 

Because of the rotational symmetry of the density operator in the spin space, the one-body correlation function in 
different spin states: 

= -^<*{(x)fc i (0)) (17) 



is instead always identically vanishing. 
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FIG. 1: Schematic view of the state of the ideal Fermi gas at zero temperature for the model considered in the Monte Carlo 
simulation, that is with N = 12 atoms and TV = 16 grid points. Each mode is a plane wave with a wavevector k = 2ns/L and 
an energy ek = h 2 k 2 /2m, where the integer s ranges from —8 to 7. The modes with \s\ < 2 are totally filled, whereas the modes 
with |s| = 3 are half-filled, the higher energy modes being empty. Ef = h 2 kp/2m is the Fermi energy. The energies are here 
in units of h 2 /rnL 2 . 




FIG. 2: Single-spin one-body correlation function g^(x) for different values of the temperature T/Tf = 1.12,0.56,0.056 
(circles, squares, diamonds). N — 12 atoms on a N — 16 points lattice. Coupling constant p^go/ksTF = —0.42. 



B. Density-density correlation functions 



Density-density correlation functions are another observable of interest. Both the single-spin density-density correla- 
tion function: 



Pa Pa Pa ax 

and the opposite-spin one: 



jffW = ^<*{(0)*{CiO*i(aO* T (°)> = T^<M*)Pl(0)>, 

i i i 4- ' I ' 4- 



(18) 



(19) 
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with pcr(x) = ^] r (x)'^ (T (x), have been calculated by Monte Carlo and plotted as a function of x respectively in 
fig|3K an d in figCJ). In fig0] we have plotted g^) (0) as a function of temperature. The magnitude of actual density 



(2) 

correlations is quantified by the difference g a J,(x) — 1. 
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FIG. 3: Single-spin (a) and opposite-spin (b) density-density correlation function (x) and gff(x) for different values of the 
temperature T/Tf = 1.12,0.56,0.056 (circles, squares, diamonds). Same system parameters as in fig|21 
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FIG. 4: Opposite spin density-density correlation function at x — 0: g^(0) as a function of the temperature T in the canonical 
ensemble. Circles: Monte Carlo results. Dotted, short dashed, long dashed, solid lines: perturbative results upto order 
respectively 0, 1, 2, and 3. Same system parameters as in figEl The vertical line is the somewhat arbitrary lower bound on the 
temperature range where perturbation theory converges rapidly. 



( 2 ) 

On one hand, the density correlations in a single spin state described by g^' show a short-range hole (Pauli hole) of 
width similar to the bump of the one-body correlation function gi 1 ) and are weakly affected by the interactions and 
by the temperature variations (figBK)- 

(2) 

On the other hand, the density correlations between opposite spins described by g^ show an interesting temperature 
dependence in the presence of attractive interactions. The lower is the temperature, the most effective are in fact the 
interactions and therefore the stronger the bunching of opposite spin particles on a given lattice site. In figOt> we 

(2) 

have plotted the spatial profile of g^,'(x) for different values of the temperature: for the lowest value of T/Tp, notice 

( 21 

not only the increase of gl, (0), but also the appearance of oscillations as a function of x. As we shall see in the next 
subsection, at this temperature a condensate of pairs is present. The oscillations then result from the contribution 
of two distinct effects: the Friedel oscillations in the correlation functions of the normal phase which follow from the 
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sharpness of the Fermi surface [l3| , and the oscillations shown by the Cooper pair wavefunction described within the 
BCS theory by the pairing function (^(^^(O)). In hgQ]we have summarized the values of '(0) as a function 

of the temperature. Notice that g?, '(0) — 1 is appreciable already at the highest temperature considered in fig0] 
which, as we shall see in the next subsection, is much higher than the critical temperature T* for the appearance of 
long-range order. 



C. First-order pair coherence function 



It is believed in statistical physics that the superfluid transition in two-component Fermi systems with attractive 
binary interactions is related to the appearance of long-range order in the so-called anomalous averages [10|. In 
symmetry breaking theories such as the BCS one, this feature corresponds to a non-vanishing value for the gap 
function defined as: 



■go(9 r {x)9i(x)), 



(20) 



which plays the role of the order parameter of the phase transition in a Ginzburg-Landau approach. In number 
conserving approaches, quantities like i|2U|) are zero. The phase transition however still appears in the long-range 
behaviour of correlation functions of the form: 
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'pair 



(x) 



1 



* T (0)*l(0)>. 



(21) 



A similar criterion was used in pll I22I l23j to determine the transition temperature. 

A simple physical interpretation of gp^ T can be provided as the first order correlation function of pairs: the operator 
*l(0)4'|(0) annihilates in fact a pair of particles in opposite spin states at the spatial position and the operator 
^{(a:) 1 I r j(:r) creates them back at x. This correlation function is therefore formally equivalent to the first order 
coherence function of a composite boson formed by a pair of fermions with opposite spins. From this point of view, 
the non- vanishing long-range limit of g^l^ix) is a signature of a quantum condensation of pairs. 





FIG. 5: Normalized pair coherence function <7p air (^) for two different temperatures T = 0.562> (left panel) and T = 0.056Tf 
(right panel) in the canonical ensemble. Circles: Monte Carlo results. Dotted, dashed, solid lines in left panel: perturbative 
results upto order respectively 0, 1 and 2 (orders 1 and 2 are undistinguishable) . Dotted, short dashed, long dashed, solid lines 
in right panel: perturbative results upto order respectively 0, 1, 2, and 3. Same system parameters as in figEl 



Monte Carlo simulations for this quantity are shown in figCO At low temperatures, <?p a i r (aO has a finite value also 
for x = L/2, i.e. at the largest distance from allowed by the finite size of the box. On the other hand, at higher 
temperatures, but still much lower than the Fermi temperature, g p ]^i r (L/2) becomes very small and the long-range 

order is destroyed. To make this cross-over more apparent, we have plotted in Fig. [5] the value of g p 1 ^ I (L/2) as a 
function of the temperature: a sudden rise of this quantity appears at low temperatures. This behavior qualitatively 
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FIG. 6: Pair coherence function g^ T (L/2) as a function of the temperature T. Circles: Monte Carlo results in the canonical 
ensemble. Solid line: BCS theory. Same system parameters as in figEl 



corresponds to the one expected for a BCS transition: although a BCS transition can not occur in one dimension in the 
thermodynamical limit because of long wavelength fluctuations destroying the long range order , it can however be 
observed in our simulations because of the finite size of the system. As the system is finite, the transition temperature 
T* is not precisely defined and the long-range order has an analytic dependance on temperature. Notice that the 

(2) 

opposite spin density-density correlation described by gl^ (0) are already important at T > T* and for T < T* they 
only get slightly reinforced. 



D. Second-order momentum space correlation function 

Another observable that has been recently proposed as a possible way of detecting the transition to a pair condensate 
state is the second-order momentum space correlation function |25| : 

G k 2) {k) = (n fcT n_ fei ) - (n m )(n. kl ), (22) 

where the operator ti% a — al^aka gives the occupation of the plane wave k with spin component a. 

As discussed in [2^, BCS theory predicts that correlations should be absent above Tbcs, thai is G^ 2 ' = 0, while the 

(2) 

transition to a condensate state should be observable as the appearance of a non- vanishing value of G k , sharply 
peaked around k = hp. 

In figd we have plotted Monte Carlo results for G? 1 as a function of k for different values of the temperature. At all 
temperatures, the quantity is indeed strongly peaked at k = fcf, and nearly vanishes at the other values. A summary 

(21 

of the temperature-dependence of G k {kp) is plotted in fig|S| At temperatures above the transition temperature T*, 
correlations are negative and increase as the temperature is lowered. The negative correlation simply follows from the 
fact that we are working in the canonical ensemble, that is at a fixed total number of particles (see the ideal Fermi 
gas result in figlBl)- As the temperature drops below T*, the correlations change sign becoming large and positive. 
This is a signature of pairing: because of the attractive interactions, the states with a filled Fermi sphere plus two 
particles in states of opposite momenta and spins are in fact energetically favoured. In this state, the fluctuations of 
the occupation numbers of the kp, | and —kp, j states are positively correlated. 
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FIG. 7: Second-order momentum space correlation function (k) for different values of the temperature. Empty circles, 
squares, diamonds: Monte Carlo results for T/Tf = 0.56,0.28,0.056. Filled circles and squares at k = ±&f: perturbative 
expansion up to order 3 for T/Tf = 0.56, 0.28 respectively. For T/Tf = 0.056 a perturbative expansion to an order higher than 
3 would be required to observe convergence. 
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FIG. 8: Second-order momentum space correlation function at k = &f: G k (&f) as a function of the temperature T. Circles: 
Monte Carlo results in the canonical ensemble. Solid line: BCS theory in the grand canonical ensemble. Dashed line: ideal gas 
in the canonical ensemble. Same system parameters as in fig|5| 



V. RESULTS OF A PERTURBATIVE EXPANSION IN g 



In this section, we explain how to calculate the pair coherence function 3p air and the density correlation function 
by means of a series expansion in powers of the coupling constant gg. The same procedure was used in fig0to obtain 

(2) 

a series expansion for the momentum-space second order correlation function G k although we do not give here the 
details of the calculation. We expect this perturbative approach performed around the ideal Fermi gas to be efficient 
mainly at T > T* that is in absence of a condensate of pairs. For T < T* we indeed found numerically that the series 
(up to order 3) is slowly convergent. Note that such a series expansion can however be shown to be convergent at 
non-zero temperature for our model system with a finite number of modes, see below. 
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A. In the canonical ensemble 



As we wish to compare to the Quantum Monte Carlo results, we have in principle to perform the perturbative 
treatment directly in the canonical ensemble with N particles. The resulting averages in the ideal Fermi gas thermal 
state are however difficult to evaluate analytically. We therefore apply the following trick to 'canonize' the grand 
canonical ensemble. We introduce an unnormalized grand canonical thermal density operator defined as 

a gc (9) = e -^(«-^) e «* (23) 

where N is the total number operator, 9 is an angle and /xo is the chemical potential of the ideal Fermi gas with an 
average number of N particles. Taking the Fourier component of <j(9) over the harmonic e amounts to projecting 
a{9) over the subspace with exactly N particles. The canonical expectation value of an operator O is therefore exactly 
given by [3l|: 



J^d9e-* dN Tr{a gc (9)Q} 
J^d9e-^Tr[a gc (9)] 



(o)n = v.: "Z - ^ 



We then expand e ^ in powers of the interaction potential V, here up to third order: 

rfl r P r T 2 



l - / dn v(n) + / dr 2 f 2 dn V(T 2 )V(n) 

Jo Jo Jo 



f dr 3 [ 3 dr 2 [ 2 d n V(T 3 )V(T 2 )V( n ) + 
Jo Jo Jo 



(25) 



where Tio is the kinetic energy operator of the gas and the imaginary time interaction picture for an operator X is 
defined as 

X(t) = e r(n a -^N) Xe -r(H -^ N)^ ^g) 

For a non-zero temperature and a finite number of grid points, the norm of the operator V(t) is finite as both V and 
T~£o — Mo-^V have a finite norm. As a consequence, the norm of the n th -order contribution to the series expansion Eq. 125|) 
can be bounded from above by A n /n\ where A is some number, and the series Eq. (|25|l is absolutely convergent |32| . 

The calculation of the numerator of Eq. (|24|l then involves the ^-dependent grand canonical partition function of the 
ideal Fermi gas and ^-dependent expectation values in the grand canonical ideal Fermi gas: 

Sq(0) = Tr[a° gc {9)} = JJ (l + e ~m 2 k 2 /(2m)~n ] e ie\ 2 (27) 



(X)o(0) = ^yTr[a g ° c (0)Al (28) 

where the square originates from the presence of two spin components. The operator X is one of the terms inside the 
square brackets of Ea. (|25|l . The expectation values can be evaluated using Wick's theorem and involve the following 
particle and hole correlation functions: 

G (x,t;#) = ${(a!,T)^ T (O)>o(0) (29) 
G (x,t;9) = $ T (a:,T)$(a))o(0). (30) 

The explicit expressions of the relevant expectation values in terms of Go and Go are given in the Appendix [BJ The 
integrals over the 'times' t\,t%,T3 and the sums over the grid points associated to each factor V{ti) are performed 
numerically. As each integral is discretized in 256 steps and there are 16 grid points in the lattice, the calculation of 
the third order correction involves the summation of about 10 10 terms for a given value of 9. 

(2) 

The perturbative results for the x — pair distribution function gW\Q) are plotted against the Monte Carlo results as 
a function of temperature in Fig. ^ for various orders of the perturbative expansion. The agreement with the second 
order expansion is perfect at high temperature, whereas the third order contribution is required to have agreement 
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at lower temperatures 33]. For a given temperature, the x dependence of glMx) predicted by the perturbative 

expansion is also in good agreement with the exact Monte Carlo results. For g^^ix), the agreement is also good, at 
high temperature in Fig. E^,, as well as at a temperature T < T* in Fig. The fact that the third order prediction is 
very close to the Quantum Monte Carlo results even when long range order is present may be fortuitous: it significantly 
differs from the second order prediction so that a calculation of the fourth order correction is required to justify the 
truncation of the series at this order. 



B. In the grand canonical ensemble 



It is actually interesting to perform also the perturbative expansion in the grand canonical ensemble: simpler analytical 
formulas can be obtained, which can be used to test existing approximate theories applicable to the grand canonical 
ensemble. The unnormalized density operator of the gas is now 

a m = e-^ n -^\ (31) 

The perturbative expansion has to be performed for a fixed value of the mean total number of particles equal to N. 
As a consequence the value of the chemical potential /i is not known in advance and has to be adjusted order by order 
in the perturbative expansion. To this end, we write 

H = |Uo + 5[i (32) 

where fi is the chemical potential of the ideal Fermi gas having on the mean a number N of particles. This amounts 
to performing the following splitting: 

H- fiN =(H - fX N) + W (33) 

where the perturbation is now 

W = V-SfiN, (34) 

both terms in W being of order gg. We shall restrict here for simplicity to a second order expansion. From Ea. i|25|) 
we get 



(/ = (O)o - Ip dnjWjnWo + f£ dr 2 £ d n (ff(#(nft + ■ ■ ■ 
l-/ fl dr 1 <W(r 1 )) + / ( fd7 a ti 2 dn(W(T 2 )W(n)) + , 



where (X) = Tr[fXgcX]/Tr[<7 gc ] stands for the expectation value in the grand canonical density operator of the interact- 
ing gas Eq. (|31|l and (X)o stands for the expectation value in the grand canonical density operator exp[— /3(Tto — [loN)] 
of the ideal Fermi gas. Expanding the inverse of the denominator in Ea. 1)35(1 and keeping terms up to second order, 
one obtains 

(O) = (O) - f drx «W(n)O»0 + f dr 2 [ ' dn {{W{n)W{n)O)) + O(g 3 ) (36) 
Jo Jo Jo 

where we have introduced the irreducible averages of products of operators A, B, C: 

({AB)) = (AB) - {A)o{B) (37) 
{(ABC)) = {ABC)o - (A) {BC)o - (B) a {AC) a - {C) (AB) a 

+2(A) a (B) a (C)o (38) 

and where we used the identity 



dniWin^o] =1 dr 2 r d Tl (W(^)) (W(ti)) . (39) 
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FIG. 9: Opposite spin density-density correlation function at x — 0: (0) as a function of the temperature T in the 
grand-canonical ensemble. Solid line: perturbative result at order 2. Dashed line: density-density RPA. Dotted line: Nozieres 
Schmitt-Rink theory. Dot-dashed line: $t^rt _ \j/\j> RPA. Thin solid line: BCS theory. The vertical line is an approximate 
lower bound on T where the second order perturbative theory is accurate; its position was determined by a comparison to 
the Quantum Monte Carlo results of Fig. 2] Same system parameters as in fig|5| The mean number of particles is fixed to 
(N) = 12. 



To calculate 5p, up to second order, we express the fact that the mean density of spin up particles is fixed in x = 0. 
As the system has translational and spin symmetry, this is equivalent to fixing the mean total number of particles. 
We therefore specialize Ea. (|36f) to the case O = t^|(0)^i-|-(0) and obtain [34| 

. 1 , Xdr 2 £dr 1 {dxfY JXl ^P2iH 2 i{P2iH 2Q P w -P 2a H 21 H w ) 

s fJ- = o.9oP + 9a — — ^ — +°K%) (40) 

ho 



n yuf i so g 

1 tfdndx^PwH! 



where dx is the spatial step of the grid, p is the total density and the following notations were introduced: 

Pij = (lp^(Xi,Ti)lpt(Xj,Tj))o (41) 
H i] = Wi x i:n)lp\(Xj,Tj))o (42) 

for integers i,j equal to 0, 1 or 2 and with the convention xq = 0,tq = 0. Note that the term of order go in <5/i 
coincides with the Hartree-Fock mean field prediction. 

In a second step, we calculate g\f(0) by taking O = ^1(0)^(0)^(0)^(0). Eliminating Sfi fr om the resulting 
expression gives: 

(p/2) 2 g\]\0) = (p/2) 2 ~g Q [ dn dx ^(Proi/io) 2 

+ 9o f dr 2 [ 2 d Tl (dx) 2 V (P 21 H 20 P W - P 20 H 21 H W ) 2 
Jo Ja Xl ,x a 

+ O(g 3 ). (43) 

In the case of a negative coupling constant g$ , both the first and second order terms are positive, leading to a spatial 
bunching of opposite spin particles, as expected for attractive interaction. 



From the comparison with the quantum Monte Carlo calculations in the canonical ensemble, we know the temperature 
range over which the second order perturbative expansion gives accurate predictions for gi, (0). For the grand 
canonical ensemble with the same mean number of particles, we expect the same conclusion to apply. We therefore 
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use a numerical integration of Ea. H43|) as a test of existing approximate theories that will be reviewed in sec I VII 
As is apparent in Fig. the density-density RPA is in very good agreement with the perturbative result, whereas 

the vpt^ — <I"]/ RPA overestimates (0) and the Nozieres-Schmitt-Rink prediction clearly underestimates it. At 
temperatures above the BCS critical temperature, the BCS theory reduces to the mean-field Hartree-Fock theory 

(2) (2) 

which gives for g^'(0) simply the ideal Fermi gas result, g^,'(0) = 1. 



VI. COMPARISON WITH APPROXIMATE THEORIES 



For the grand canonical ensemble, several approximate many-body theories exist which can be used to obtain pre- 
dictions for the correlation functions of the interacting Fermi ga s. In the next subsections, some among the most 
famous ones are discussed, namely the mean- field BCS theory jiaUiLtwo versions of the random-phase approxima- 
tion (RPA) [HEl, and tne Nozieres-Schmitt Rink (NSR) theory @] developed to study the BCS-BEC crossover 
in strongly interacting gases. A quantitative comparison with the prediction of a grand canonical version of the 
perturbative expansion of section^will be performed. In order for the comparison to be meaningful, the many-body 
theories under investigation have been specialized to the specific case of the discrete lattice Hamiltonian Q with 
exacly the same discretization parameters as used in the previous sections. 



A. BCS theory 

In the BCS theory, the equilibrium density matrix is determined in a self-consistent way from the mean-field quadratic 
Hamiltonian of the grand canonical ensemble at a chemical potential fi: 



^bcs = E - m)4<A* + g J2 dx P-<rH(x)**(x) - E (*t(*)*i(s)A* + (a:)A) (44) 



'h 2 k 2 

A 

ka 

where the mean density p a in a given spin component and the gap function A are defined as usual as: 

Pa = (*Ux)9„(x)) (45) 
A = - 9 o(*tW*lW)- (46) 
The quadratic Hamiltonian (|44|l is easily diagonalized by a Bogoliubov transformation in the plane wave basis: 

Wbcs=E^4A^ ( 4? ) 

kcr 

where the c, & operators satisfy Fermi anticommutation rules and are related to the Fermi field operators by: 

^ = T^E"^ £ki+vke ikx St fe>T (48) 



*T (*) = ^= E u * e ' lkX ^ ~ "* ^ £ -H (49) 



where the positive coefficients U}~, Vk of the Bogoliubov transformation are defined by: 



< = 1 - V * = i[ 1+ E^ > (50) 



the quasi-particle energies are given by: 



h 2 k 2 ' 2 



^ = v A2 + (^r-A) (5i) 
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and the chemical potential is shifted as p, = ji — gap a so as to take into account the mean- field energy. The self- 
consistency equation for the gap A is: 

where the f k are the quasi-particle occupation numbers fk = (e Ek ' kBT + l) -1 . For high temperature T > Tbcs, 
Eci. (|52[1 has no solution, so that the system is in the normal phase A = and the BCS theory reduces to a Hartree- 
Fock theory. At low temperature T < TbcSj the gap equation is solved for a non- vanishing value of A. This value 
grows as the temperature decreases. 



1. Calculation of correlation functions within the BCS theory 



The expansion of the field operator l|49l) in terms of quasi-particle creation and destruction operators can be used 
to obtain a prediction for the correlation functions. For instance, the BCS prediction for the one-body correlation 

function <jw(#) is given by: 



9*K X ) = — ^'—j—[\ u k\ 2 h + \vk\ 2 {l-fk)\- (53) 
Pa k L 



(2) 

From Wick's theorem, the single-spin density-density correlation function g&J (x) is 

g%{x) = l-\g£}{x)\\ (54) 

Both quantities are affected in a weak way by the attractive interactions and eventually by the appearance of a 
non- vanishing gap A. 

A richer physics can be found in the opposite spin density-density correlation function gi/(x). For this quantity, the 
BCS theory predicts: 



g\ z >(x) = l + \A(x)\ z , (55) 



ii 2 hx) = i + ^u(x)i 2 

where the anomalous correlation function A{x) is defined as: 

A(x) = (*iOr)* T (0)) = Y, ^—ukv k (l - 2/*). (56) 



As A = — goA(O), 0ir (O) has the simple expression: 



0f>) = i 



A 2 



90 Pa 



(57) 



For T > Tbcs 5 this quantity is identically 1, which means that the BCS theory does not predict any correlation 
between the densities in opposite spin states. These appear only for T < Tbcs as a consequence of the non-vanishing 
BCS gap. As one can see in figEI this result is in qualitative disagreement with the perturbative expansion which 
gives a significant degree of correlation also for T ~ TV ^S> Tbcs ■ 

The BCS prediction for the first-order pair coherence function g^ir( x ) i s: 

9%{x) = -^\^\ 2 + &Xx? (58) 
P^Pi% 

and is characterized by a short-ranged bump of spatial size of the order olip = 1/kp, and a non- vanishing long-range 
limit. As one can see in figEl the long distance behaviour of gpaL predicted by the BCS theory is in qualitative 
agreement with the Monte Carlo predictions. 
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B. Random Phase Approximation 



1. Fluctuation-dissipation theorem 



A simple way of including the fluctuations around the mean-field is to compute a response function within the 
mean-field theory and then invoke the fluctuation-dissipation theory to obtain the corresponding correlation function. 
In this subsection, we shall give a short review of the main results of linear response theory that are required to 
obtain the correlation functions of our interacting Fermi gas. A complete discussion of linear response theory and 
fluctuation-dissipation theorem can be found in |ld l26| . 

Let A and B be two operators of a system characterized by a time-independent Hamiltonian TL. For notational 
simplicity, we assume that at equilibrium (A) eq = (B) oq — 0. A weak perturbation of the form: 

H P ert = e{t)A^ + e*{t)A (59) 

is applied to the system and its effect on the observable B recorded. At linear regime, this is summarized by the 
linear response functions: 



(B)(t)= df XBA (t-t'y{t') + x B Ai{t-t'W) ■ (so) 

J — oo 

The linear susceptibilities x have the simple expression in terms of commutators: 

XBA(t) = ^Tr{[B(t),A]p eq }Q(t), (61) 

where p oq = ~g e ~ l3 ' H is the thermal equilibrium density matrix at (3 — 1/ksT, Z — Tr[e~^ w ] is the partition function 
and B(t) — e int / h Be~ lUt l h . As the Hamiltonian of the unperturbed system does not depend on time, the Fourier 
transform of XBA(t) is the frequency-dependent response function to a harmonic perturbation of frequency w: 



X baH= / dte^XBA^e-^, (62) 

J — OO 

where 77 — > + [3|| The correlation function Sba{^) is denned as: 

p OO pOC 

S BA {oj)= dte wt S BA (t)= dte wt (B{t)A{Q)). (63) 



If the condition: 

Tr[V En BV Em AV En ]e M, (64) 

holds for all the eigenenergies E n and E m of the Hamiltonian TL where the projector Ve projects onto the eigenspace of 
energy E, then the fluctuation-dissipation theorem holds in its most common form (Callen-Welton theorem) relating 
the imaginary part of the response function Im^B^w)] to the correlation function Sba{u): 

lm[ XBA (cu)] = -±.§ BA (u)(l - e~^). (65) 

It is easy to verify that the condition (|64|l is verified if A* = B or, more generally, if A' = SB S, S being an arbitrary 
unitary operator such that S 2 = 1. 

The fluctuation-dissipation theorem H65[) implies that the correlation function Sba (w) is fixed by the knowledge of 
Im[xsA(w)] modulo a delta distribution in uj = 0: 

Sba(u) = -2 ^^-/^ + 2nCBA ( 66 ) 
The constant Cba can be written as follows: 

Cba = -k B T[ X %T u - lim Xba M] (67) 
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in terms of the thermodynamic (isothermal) susceptibility Xba™- As usnal m thermodynamics, this is defined as the 
response on B when the system is at thermal equilibrium in the presence of a weak and time-independent perturbation 
Hport =e*A + eA^. 

SB thc rm = £ * x thcrm + e x ther m (6g) 

resulting from the expansion to first order in e of: 

X nthcrm _ ±L L" ° J ((?q\ 

Notice that while the definition of XbI"" involves some implicit coupling to a thermal reservoir at temperature T, 
Xba{<^) is defined for an isolated system evolving under the Hamiltonian H.. For this reason, the thermodynamical 
susceptibility Xsl™ an d the static limit xba{u ~ * 0) in general are not equal p^ . 

From the microscopic expression of Cba in terms of the eigenstates of 7i of energy E n : 

n 

one concludes that Cba is not vanishing in the presence of degeneracies between the eigenenergies of TL or when the 
diagonal matrix elements (n| A \n) and (n\ B \n) are not vanishing. 

The possibility of having a term in 5{uj) in the correlation function Sba is often neglected in statistical mechanics 
textbooks, e.g. |lflj . Although this is generally correct in the thermodynamical limit, it may lead to incorrect results 
for the correlation functions of finite systems. Examples of this issue are discussed in the next section and in the 
Appendix [CJ 



2. Density- density RPA 



A prediction for the density-density correlation function of the Fermi gas in the grand canonical ensemble can be 
obtained by applying the general results of the previous subsection to the operator p a (x) giving the particle density in 
the spin state a at position x. An approximate prediction for the density-density response functions can be obtained 
by linearizing the equations of the mean-field theory discussed around the thermal equilibrium state. For historical 
reasons, this approximation scheme is usually called random phase approximation (RPA) 0,0. For the sake of 
simplicity, we shall limit ourselves to the case T > TbcSi regime in which the vanishing of the anomalous averages 
considerably improves the physical transparency of the formulas. 

As the system is spatially homogeneous with the same density p a in each spin-component, the different Fourier 
components of the spatial density 

^ = "Tr E dxe ~ lkX M*) ^ P°) = -7g E dx e ~* kX (^U^a(x) - Pa ) (71) 

are decoupled. Taking A = dpta and B = Sp k , a ,, the frequency-dependent susceptibility matrix has the form: 

Xcrcr' (fc, k'; lo) = Xaa> (k, uj) 8 k ,k' ■ (72) 

Because of the symmetry in the spin space, XTT = XU an d XT J. = Xlh so that the eigenvectors of the susceptibility 
matrix are the symmetric and antisymmetric linear combination of the two spin states. The corresponding eigenvalues 
are: 

X±{k,u) = Xu{k,uj) ± xn(fc,w). (73) 
Conversely, the susceptibility matrix x<t<j' in the a =\\, basis is written as a function of the x± a s: 



1 



x+ + x- x+ - X- 



2 V x+ ~X- X++X- 



(74) 
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The RPA susceptibility of an interacting gas can be calculated from the Hartree-Fock equation of motion [27|. and 
has the simple expression: 

X±(k,uj ) = — -r-77 — s, (75) 

lT5oXo(fc,^) 

in terms of the susceptibility xo of a non-interacting, one component Fermi gas at the same temperature and chemical 
potential: 

L nui + c q — c q+ k + i 0+ 

£ q = h 2 q 2 /2m are the energies of the single-particle states and f q = (1 + exp[f3(£ q — /i)]) _1 the corresponding Fermi 
occupation factors. Notice that the quantity inside the sum vanishes for the states q such that £ q = £ q +k- 

An expression analogous to Ij75(l : 

therm/ u\ 
therm/, \ _ AO VV 

relates the thermodynamic susceptibilities x± eTm (k) of the interacting gas to the thermodynamic susceptibility of the 
non-interacting one-component gas: 

! e g(£ g -£ g+fc )_ 1 -rc/c 
U &q /= Cg+fc 
(78) 
/3 if £g = £ ?+fc 

By applying the fluctuation-dissipation theorem in its form <|6(j|) to the operators B — 6pu a and A — Spka' , one can 
write the correlation function 

/oo 
dte wt {5p{ a {t)5p ka ,{Q)) (79) 
-OO 

in terms of the imaginary part of Xaa' (&> ^) an d the thermodynamic susceptibility X^ 1 J r c / rm (A;): 

= ~2fi Im i ^'_ ( ^ )] ~ 2^fc B T[ x t CT h CT 7 m (fc) - lkw(k )W )] (80) 

The condition ffity is here satisfied since A and B are connected by £?t = SAS with 5 respectively equal to the identity, 
if a = a', or the spin-inversion operator S exchanging the spin components f , J, of all the particles, if a = — a' . 

Finally, the RPA prediction for the desired real-space, one-time density-density correlation function can be found by 
inverse Fourier transform of S aa i(k,uj): 

(p a ( x )p al (Q))=p aPa , + - J —J2e- ikx S„„,(k,u;). (81) 

Corresponding predictions for the opposite spin density-density correlation function at x = as a function of the 
temperature are plotted in figGD Notice the excellent agreement of the RPA prediction with the one of the perturbative 
expansion in go discussed in sec0 In our finite system, the agreement strongly relies on the correct inclusion of the 
S(ui) term in (|80|) . An explicit calculation of this issue for the non-interacting case is presented in the Appendix [U] 



3. $t^t.$q, upa 



In the previous subsection, we have obtained a prediction for the density-density correlation function g a ^, (x) of an 
interacting Fermi gas by using the RPA density-density susceptibility and then invoking the fluctuation-dissipation 
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theorem. In the present subsection, a similar approach is used to obtain the pair coherence function g^ ir (x) at 
temperatures higher than the BCS critical temperature Tbcs m the grand canonical ensemble. 

Consider the pair of operators: 

B = §\{x)i>\{x) (82) 
A = $ T (0)#;(0). (83) 

At thermal equilibrium both of them have a vanishing expectation value. The correlation function: 

(BA) = ($\(x) %\(x) * T (0) * ; (0)) = P1PI g£L(x) (84) 
can be evaluated from the susceptibility xba- 

As in the previous subsection, we introduce the spatial Fourier components as: 

A k = ^]r^e-^ T (x)^(x) (85) 



1 

7t 



B k = ^=^dxe 4fca; *{(x)*{(x); (86) 



note the sign difference in the phase factors of (|85|l and (|86[) . Thanks to the spatial homogeneity of the system, 
the susceptibility is diagonal in fc-space. From the Hartree-Fock-Bogoliubov equation of motion for the anomalous 
averages |27j . one can obtain the following simple expression for the RPA susceptibility: 



where xo is defined as: 



X{k,Uj) = - — ; r, (87) 

1 - 9oXo{k,u) 



L ^ lj + £ q + £ q+k + g a p — 2p + i0+ 



and describes the ideal gas response. As previously, £ q are the energies of the single-particle states, f q the Fermi 
occupation factors, and p = p^+pi the total particle density summed over both spin states. Notice how x(k = 0, lo = 0) 
diverges when x.o(k = 0, u> = 0) tends to l/go- This is the signature of the approaching of the BCS transition: the 
standard equation ^(J for the BCS critical temperature is in fact recovered if one imposes: 

i = . 9 o xo(k = o, u, = o) = -|| £ 1 r^rr .. - ( 89 ) 



'I 



q + 2.9o« - P 



As the chemical potential is a variable that can be continuously varied, all degeneracies between states with different 
particle number are accidental and occur only for discrete values of p. As A and B have vanishing diagonal elements, 
the correction term in 8(u>) vanishes for all other values of p. As <7pa; r has a continuous dependance on p, there is no 
need for calculating x thcim (k)- We therefore have: 

S BA {k, u) = -2h - _ — (90) 

and 

~fjE e-^S BA (k,u)- (91) 

PWl L J-OO I™ 

The condition (|64|l is here satisfied as B\ — A k . As g^ir( x = 0) coincides with g?, '(x = 0), we have included in figEl 
also the prediction of the present vpttjrt _ vj/ij/ RPA approach. The agreement with the perturbative expansion is less 
good than in the case of the density-density RPA approach. 



19 



C. Nozieres-Schmitt Rink approach 



In @, a non-perturbativc calculation is performed for the grand potential of a two component Fermi gas with attractive 
interactions by a resummation of a certain class of diagrams. 

Starting from Eq. (20) of @ which gives the grand potential n in terms of an integral in the complex plane, we can 
deform the integration contour and apply the residues formula to obtain the following expression in terms of a sum 
for the case of a contact interaction potential 36] : 

00*, T; g ) = n (M, T) + k B T £ log[l - x (q, u v ;go)] ~ E t 1 ~ ^ - ( 92 ) 

q,u v ki,k 2 

where w v = TA/kv j ' f3 with v integer ranging from — oo to +oo. The function \ is defined as: 

I \ 9° 1 - fk - fk+q 

X(q,u;g ) = -— > - — * — , (93) 

k k k+q ~ P ~ u 

f q being the occupation number f q = {l + exp[f3(£ q — //)]} _1 and FIq being the grand potential for the ideal two- 
component Fermi gas. 

This theory requires a self-consistent determination of /i from N = —d^tt, which we perform numerically. It gives 
access to g^ (0) thanks to the Hellmann-Feynman theorem [2^ : 



The results are plotted in Fig|5J The poor agreement with the perturbative expansion can be explained as follows. 
Let us expand Eq. l|92l) upto second order in go at a given fi. If one then replaces /x by its value in the NSR theory 
for the density under consideration, one gets from Ea. (|94|l a prediction for g^ (0) upto first order that the can be 
compared to the exact expansion Ea. (|43|) : 

i&(0) - 9^(0) = P ^>f^ 1 + 0{gl) (95) 

Here po(/ i ) is the total density of the ideal two-component Fermi gas for the chemical potential /i. A first order 
expansion of £1 is enough to obtain the chemical potential Unsr = Mo + 5op/2 + 0{go) i n the NSR theory and to 

(2) 

conclude that g NSR differs from the exact value by a term of the order of g which has the same sign as g . This just 
because not all the second order diagrams for n have been included in the resummation procedure |37j . 



VII. CONCLUSIONS 



In the present paper, we have presented the result of extensive Quantum Monte Carlo simulations for the static 
correlation functions of a one-dimensional lattice model of attractively interacting two component fermions. The 
numerical results obtained by QMC have been compared to existing approximate theories. Excellent agreement with 
the predictions of a perturbative expansion in the interaction constant has been found, as well as with the ones of the 
random phase approximation. 

Although long-range order is destroyed by phase fluctuations in one-dimensional systems in the thermodynamical 
limit, the finite size of the system under consideration still allows for the identification of a crossover to a condensed 
state at the temperature T* at which the first order coherence length of the pairs becomes larger than the system 
size. 

We have found that a significant degree of opposite spin density-density correlations already exists at temperatures 
well above T* and is only slightly enhanced as the temperature goes below T*. This means that a measurement of 
the density-density correlation function gl/ can not provide an unambiguous signature of the onset of a condensed 
state in the gas. On the other hand, this could be provided by a measurement of the second-order momentum space 
correlation function as suggested in pjj . or, even more directly, of the long-range behaviour of the first-order pair 
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coherence function g^ ir (x). A non- vanishing limit of gQ T (x) for large x corresponds in fact to the presence of a finite 
condensate fraction in both weak- (BCS) and strong- (BEC) interaction regimes. A possible experimental scheme to 
measure in atomic Fermi systems by means of matter-wave interferometric techniques will be the subject of a 
forthcoming publication. 
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APPENDIX A: THE MONTE CARLO SAMPLING ALGORITHM 

In the present appendix, we describe the numerical algorithm used for the numerical simulations. A Monte Carlo 
technique has been used to sample the probability distribution of the initial wavefunctions and the elementary 
noise terms (x) . 

As the effective contributions of the different realizations to the observables involve the trace of the Hartree-Fock 
ansatz a, i.e. the scalar product between the two A^-body Nartree-Fock states, they can have enormously different 
values, so a direct draw of the random variables (f>f^ and Q (x) would be poorly efficient. An importance sampling 
scheme [2^ has therefore been implemented, using the value of the modulus of the trace at the end r = (3 of the 
imaginary-time evolution as the a priori probability distribution function P : 

^o[{0f ) ,ej Q) (^)}] = |Tr[cr]| = |{^ 2) (/3)...^(/?)|0i 1) (/?)...^ ) (/3))|. (Al) 

In this way, the contributions of the different realizations to the trace have the same absolute value, although their 
phases are still random. 

In order to sample Pq, a Metropolis scheme |30) has been implemented: at each step a random move is proposed for 
both the wavefunctions </> l -°' ) and the elementary noises ^ a) (x). For the first one, a rotation lZ n g in the one-body 
Hilbert space is chosen with random rotation axis n and angle 9, and then is applied to all the orbitals (fif*^: 

ft (0) = n n9 (A2) 

This operation rotates the hyperplane spanned by the set of orthonormal orbitals {</>^} describing the initial Hartree- 
Fock state. For what concerns the elementary noises ^ a \x), each of then is displaced to a new position Q a '(x) as 
follows: 

^ a \x) = yw^w + ^f'w, (A3) 



b " {x) being independent, zero-mean, complex Gaussian variables such that bj (x) 2 = 0, \bj (x)\ 2 = 1. This kind 
of random process is such that the resulting distribution of the £'s is indeed a Gaussian with the required width. 
The parameter 77 as well as the probability distribution for the random rotation angle 9 are free parameters which 
can be tuned to optimise the efficiency of the simulation. Denoting with Pq" 1 ^ and Pp fin ' the value of the a priori 
probability for the configurations respectively before and after the proposed move, this is accepted with a probability 
p = min[l, Pg fin ^/Pg ln ^]. As all configurations can be attained by the random motion and detailed balance is verified, 
the stationary probability distribution of the stochastic process is indeed the desired one Pq . If a large enough number 
of moves is performed between successive realizations, these can be considered to be statistically independent. 
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APPENDIX B: EXPECTATION VALUES FOR THE CALCULATION IN THE CANONICAL ENSEMBLE 

The calculation of g^ (x a ) involves the following expectation values: 

ln{0) = (Pt(Xn,r n )Pl( X n,T n ) ■ ■ ■ P\ (x , T )/3j. {x a , T a )} (0), (Bl) 

where p a (x,r) = ^(x, t)Vv(x, t) and x = 0, r = r a = 0. Introducing the notations for the particle and hole 
correlation functions 

Py = (^ T t (x l ,r 4 )V- T (x,,T,)) o (0) (B2) 

Hy = (V' T (x 1 ,t 4 )V-{(x,,t,)) o (0) (B3) 

where i,j are integers from to n or are equal to a, we find 

Io = i& (B4) 
h = (P 2 o + PiaH la )(a -> 0) (B5) 

h = + P2lH 21 P a0 +PlaH la P a0 +P2aH 21 H la 

-P2iH 2a P la + P 2a H 2a P OQ )(a^0) (B6) 
h = {P32H 32 P$ a + P 32 H 32 P la H la - P 32 H 31 P 2a H la 

— P 3 2H 3a P 2a PoQ + P 32 H 3a P 2 iPi a — ^32^31 -P2I-P0O 
+-P2a-f^21-ffla-P00 + ^o^a-PoO ~~ ^31-f^3o-PlofoO 

+ p3 aH 3 \H\ a P — P3aH 3 iP 2 iH 2a + P 3 lH 32 H 2 iP 00 

— P3lH 32 H 2a Pi a + P 3a H 3a PQ + P 3a H3aP2lH 2 i 

+P 31 H 31 P 2a H 2a + P 3 iH 3l P^ + P 4 + P la H la P^ 
+P 2 iH 2 iPq — H 2a P 2 iPi a Poo + P 3a H 32 H 2 iHi a 

+P 3a H 32 H 2a P 00 - P 31 H 3a P 2a H 21 )(a -> 0) (B7) 

The property that each of these expressions is a product of two similar factors, the second one deduced from the first 
one by the replacement a — > 0, originates from the fact that the ideal Fermi gas in the grand canonical ensemble 
consists of two equivalent and independent spin components. 

Similarly the calculation of g^ ir {x a ) involves the expectation values 

)...${x a )4\{x a )$ l (p)fa(p)). (B8) 

Using the previous notations we then obtain 

Jo = Plo (B9) 
Ji = (PooPao + P W H la f (B10) 

'h = {P2oH2lHla + H 2a P 20 PoO + Po()PaO 

+P 10 H la P 00 + P a0 P 2 iH 21 - P 21 H 2a P w ) 2 (Bll) 
J3 = {P3iH 3 \P a oPoo + P OQ P a o — P 3 iH 3a P\oPQQ + P m P\oH\ a 

+ P 2 oH 2(1 Pq Q + P2C>Pl2lHi a Poo — P 2 lH 2a PloPoO 

+P 2 iH 2 iP 00 P a0 + P 30 H 3 iHi a P o + P 3 $H 32 H 2a P m 
+P 32 H 32 P a0 P o — P3 2 H 3a P 2 nPoo + P3iH32H 2 iP a o 

— P 32 H 3 \P 2 iP a a + P 32 H 32 PioHi a — P 32 H 3 iP 2 oHi a 

+P 32 H 3a P 2 iPio + P 3 oH 32 H 2 iHi a — P 3 oH 3 iH 2a P 2 i 

+ P3oH 3a P 2 lH 2 i — P31P132H2C1PW + ^31^31^20-^20 

—P3iH 3a H 2 iP 20 + P 3 oH 3(1 Pq ) 2 . (B12) 

The fact that these expressions are squares is again due to the existence of two independent and equivalent spin 
components. 
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APPENDIX C: DENSITY FLUCTUATIONS OF AN IDEAL FERMI GAS: PHYSICAL MEANING OF 

THE 5(u) TERM 



The density-density correlation function for a one-component non-interacting Fermi gas in the grand canonical en- 
semble at inverse temperature (3 = 1/ksT and chemical potential fj, can be calculated from Wick's theorem: 

G&(x) = (p{x) p(0)) -P 2 4E e<te /«( J ~ ( C1 ) 

qk 

where p(x) is the operator giving the density at point x. A calculation based on the fluctuation-dissipation theorem 
neglecting the term 2itCba S(uj) term in (jfifij) would give: 

qk 



where the sum has to be performed over the pair of states such that £ q 7^ £ q +k- Using the relation: 

e /3(£,-M) 5 ( C3 ) 



1 fq _ p P(S q -n) 



h 

one can see that expression l|C2(l does not coincide with (|C1J| because of the missing contribution of the pairs such 
that £ q — £ q +k- 

Inclusion of the term proportional to S(u)) in l|6tj[) fixes the problem, since it exactly provides the missing contribution: 



AG™(x) = -]T e ik *f q (l f g+k ), (C4) 

qk 

where the sum nas to be performed over the pairs of states such that £ q — £ q +k- The physical meaning of the 
k = term which contains the contribution of the diagonal matrix elements of the perturbation is transparent: it 
keeps track of the total particle number fluctuations of the grand canonical ensemble. In our spatially homogeneous 
system, degenerate pairs of states for k = — 2q are also present. 

It is apparent from (|C4|I that the contribution of the AG' 2 ) correction term tends to zero in the thermodynamical 
limit L — > 00. 
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